home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_8 / a8_1.m next >
Encoding:
Text File  |  1994-06-05  |  2.8 KB  |  117 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 8.1 (Golden Search for a Minimum).
  9. % Section    8.1,    Minimization of a Function, Page 413
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program finds the minimum of a function f(x).
  13.  
  14. % The method is the golden ratio search.
  15.  
  16. % It is assumed that f(x) is unimodal over [a,b].
  17.  
  18. % The function  f(x)  is to be placed in the M-file  f.m
  19. % function z = f(x)
  20. % z = x.^2 - sin(x);
  21.  
  22. delete f.m
  23. diary f.m; disp('function z = f(x)');...
  24.            disp('z = x.^2 - sin(x);');...
  25. diary off;
  26.  
  27. f(0); % Remark. f.m and golden.m are used for Algorithm 8.1
  28.  
  29. pause % Press any key to continue.
  30.  
  31. clc;
  32.  
  33. % Place the endpoints interval [a,b] in  a  and  b.
  34.  
  35. a  = 0;
  36.  
  37. b  = 1;
  38.  
  39. pause    % Press any key to graph the function.
  40.  
  41. clc; clg;
  42. hs = (b-a)/150;
  43. Xs = a:hs:b;
  44. Ys = f(Xs);
  45. axis([-0.05 1.05 -0.25 0.16]);...
  46. plot(Xs,Ys,'g');...
  47. hold on;...
  48. plot([-0.05 1.05],[0 0],'b',[0 0],[-0.25 0.16],'b');...
  49. xlabel('x');...
  50. ylabel('y');...
  51. title('To search for the minimum of y = f(x).');...
  52. grid;...
  53. axis;...
  54. hold off;...
  55. shg; pause    % Press any key to continue.
  56.  
  57. clc;
  58. % Place the endpoints interval [a,b] in  a  and  b.
  59.  
  60. % Place the tolerance for the abscissas in  delta.
  61.  
  62. % Place the tolerance for the ordinates in  epsilon.
  63.  
  64. a  = 0;
  65.  
  66. b  = 1;
  67.  
  68. delta = 1e-5;
  69.  
  70. epsilon = 1E-7;
  71.  
  72. [p,yp,dp,dy,A,B,C,D] = golden('f',a,b,delta,epsilon);
  73.  
  74. pause    % Press any key to find the minimum of f(x).
  75.  
  76.  
  77. clc; clg;
  78. hs = (b-a)/150;
  79. Xs = a:hs:b;
  80. Ys = f(Xs);
  81. axis([a b -0.25 0.175]);...
  82. X1 = A(1:7); Y1 = f(X1);...
  83. X2 = B(1:7); Y2 = f(X2);...
  84. X3 = C(1:7); Y3 = f(X3);...
  85. X4 = D(1:7); Y4 = f(X4);...
  86. plot(Xs,Ys,'g',X1,Y1,'or',X2,Y2,'or',X3,Y3,'or',X4,Y4,'or');...
  87. hold on;...
  88. plot([a b],[0 0],'b',[0 0],[-0.25 0.175],'b');...
  89. xlabel('x');...
  90. ylabel('y');...
  91. title('The golden search for the minimum of y = f(x).');...
  92. grid;...
  93. axis;...
  94. hold off;...
  95. shg; pause    % Press any key to continue.
  96.  
  97. Mx1 = 'The results for the golden ratio search.';
  98. Mx2 = 'The solution is [p  f(p)];';
  99. Mx3 = 'The  abscissa  p ';
  100. Mx4 = 'error bound for p';
  101. Mx5 = 'The minimum value f(p)';
  102. Mx6 = 'error bound for f(p)';
  103. clc,echo off,diary output,...
  104. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(''),...
  105. disp(Mx3),disp(p),disp(Mx4),disp(['± ',num2str(abs(dp))]),disp(''),...
  106. disp(Mx5),disp(yp),disp(Mx6),disp(['± ',num2str(abs(dy))]),diary off,echo on
  107.  
  108. pause % Press any key to see all the iterations.
  109.  
  110. Z = [A;C;D;B]';
  111. Mx7 = 'The complete list of iterations:';
  112. Mx8 = '     a(k)      d(x)      d(k)       b(k)';
  113. clc,echo off,diary output,format short,...
  114. disp(''),disp(Mx7),disp(''),disp(Mx8),disp(Z),diary off,echo on
  115.  
  116.  
  117.